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The dynamics of a system formed by a finite number of globally coupled bistable oscillators 
and driven by external forces is studied focusing on a global variable defined as the arithmetic mean 
of each oscillator variable. Several models based on truncation schemes of a hierarchy of stochastic 
equations for a set of fiuctuating cumulant variables are presented. This hierarchy is derived using 
Ito stochastic calculus, and the noise terms in it are treated using an asymptotic approximation 
valid for large A'^. In addition, a simplified one- variable model based on an effective potential is also 
considered. These models are tested in the framework of the phenomenon of stochastic resonance. 
In turn, they are used to explain in simple terms the very large gains recently observed in these 
finite systems. 



PACS numbers: 05.40.-a,05.45.Xt 



I. INTRODUCTION 



Noise induced phenomena in nonlinear systems have 
attracted a great deal of attention in a variety of contexts 
in physics, chemistry and the hfe sciences. An impor- 
tant example is the phenomenon of stochastic resonance 
(SR) [T], in which the response of the system (output) 
to external driving (input) is amplified and optimized for 
certain values of the noise parameters. More specifically, 
the non-monotonic behavior of the output signal-to- noise 
ratio (SNR) with the strength of the noise is a widely 
accepted signature of SR. In addition, a dimensionless 
quantity known as the SR gain is usually defined as the 
ratio of the output SNR over the input SNR. 

Very recently, very large SR gains have been reported 
for systems formed by a finite number A'^ of globally 
coupled bistable oscillators 0, Here the term global 
coupling is used to indicate that each oscillator inter- 
acts with all other oscillators. These systems were used 
years ago by Kometani and Shimizu ^4] as an empirical 
model to describe muscle contraction. Later on, Desai 
and Zwanzig [sj gave a more detailed statistical mechan- 
ical description, finding an order-disorder transition for 
a variable defined as the expectation value of the posi- 
tion of one oscillator. This variable is used to study the 
global behavior of the coupled bistable system. Desai and 
Zwanzig focused on systems with infinitely large sizes, 
N ^ OO, investigating the system dynamics by analyz- 
ing the time evolution of the central cumulant moments 
of one-oscillator's distribution. In addition, a Gaussian 
approximation was proposed in order to close the cumu- 
lant moment hierarchy and obtain analytical expressions. 
A similar approach is currently used as a mean field ap- 
proximation in the investigation of various noise-induced 
phenomena such as noise-induced phase transitions 
Recently, in order to study the effect of fluctuations due 
to the finite size of the system, Pikovsky et al. ex- 
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tended this approach by replacing the expectation val- 
ues of one-oscillator properties (•) by arithmetic means 

over all oscillators A^^^X]i!=i(')- ^ Gaussian approxi- 
mation, including noisy terms, was derived and used to 
illustrate the phenomenon of system size resonance, in 
which the SR quantifiers display a non-monotonic be- 
havior as a function of N. In this paper, the work by 
Pikovsky et al. is extended to higher orders in the fluc- 
tuating cumulant dynamics. The Gaussian approxima- 
tion is re-derived using a rigorous formalism based on Ito 
stochastic calculus and compared with other approxima- 
tions. 

One important goal of this paper is to explain the very 
large gain values observed in globally coupled bistable 
systems |2|, Q , especially when compared with those ob- 
served in uncoupled or isolated bistable systems. To 
that effect, it is desirable to derive a simplified theory in 
which the number of degrees of freedom is much smaller 
than the number of coupled oscillators, thus being more 
amenable to analytical treatment or qualitative interpre- 
tation. In this regard, the Gaussian approximation is a 
practical alternative, though in principle not fully sat- 
isfactory, because it is not based on a small parameter 
expansion but on an uncontrolled assumption (the ne- 
glect of cumulants higher than the second) that is known 
to be not accurate even in the limit of an infinite system 

i. 

In this paper, this approximation is presented, as well 
as other simplified models with a reduced number of de- 
grees of freedom which are able to mimic the most impor- 
tant features of the system dynamics with a finite size. 
These simplified models represent different approxima- 
tion schemes and might be regarded as an expansion or 
generalization of the work by Desai and Zwanzig Q and 
Pikovsky et al. [Sj]. 

The paper is organized as follows. In the next section, 
the model system and the SR quantifiers are defined. The 
simplified models are presented in Sec. lIIII These models 
are compared to the original model system by means of 
computer simulations in Sec.|lVl Finally, SecFVlprovides 
a short summary and conclusions. 
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II. MODEL AND DEFINITIONS 



r2 = 27r/T being the driving frequency, and 



Let us consider a set of N interacting bistable oscilla- 
tors, each one of them characterized by a single degree of 
freedom Xi {i — 1, . . . ,N), whose dynamics is governed 
by the Langevin equations 0, @ 

e ^ 

= - xl + — Y^i^J - + ^rit) + F{t), (1) 

i=i 

where ^i(t) is a Gaussian white noise with zero average 
and 



m'Kj{s))=2DS,,6{t-s), 



(2) 



is a coupling parameter defining the strength of the 
interaction between oscillators, and F{t) is an external 
driving force of period T. 

To characterize the system as a whole we define the 
collective or global variable S{t) as 



1 ^ 



(3) 



The stochastic resonance quantifiers for this variable are 
defined in the usual way. The one-time correlation func- 
tion, 



C{r) = ^ l\t(Sit + r)S{t))._ 



(4) 



can be written as the sum of two contributions: a coher- 
ent part, 

Ccoh{r) = ^ [ dt{S{t + T))^{S{t))^, (5) 
^ Jo 

which is periodic with period T, and an incoherent part 

CincohM =C(t) -Ceoh(r), (6) 

which decays to zero for large values of r and reflects the 
correlation of the process S{t) at different times due to 
fluctuations. In the expressions above, the notation (■ • ■) 
indicates an average over the noise realizations and the 
subscript "cxd" indicates the long time limit of the noise 
average, i. e., its value after waiting for t long enough that 
transients have died out. The SNR of a random signal 
measures the signal strength relative to its background 
noise. More specifically, we calculate the output SNR as 



Rn 



~Qi 



where 



2 

Qu = 7f dr Ccoh(T) cos(r2r). 



(7) 



(8) 



Qi = - 



dr Cincoh(T) cos(fiT). 



(9) 



Note that the quantity Qu is proportional to the so-called 
spectral amplification, which is another widely used SR 
quantifier. 

As the size of the system N is increased while keep- 
ing the interaction parameter 6 constant, the collective 
variable S{t) becomes less noisy, becoming completely 
deterministic in the limit N —>■ oo. As a result, i?out 
diverges in that limit. This is a consequence of the aver- 
aging process implicit in the definition ([3]). 

The SR gain is defined as 



G 



Ro 



Rh 



(10) 



where Rm is the SNR of the collective input signal 

example, for a periodic rect- 
angular driving force of amplitude A, the input SNR is 
given by i?i„ = iA^N/inD). The SR gain HUD is a di- 
mensionless quantity that measures the amplification of 
the system response with respect to the collective input 
signal. The input SNR Rin diverges linearly with N in 
the limit N — > oo, so that the SR gain remains finite. 

Since in a system with coupled linear oscillators the 
SNR of the collective process equals the SNR of the col- 



lective input signal, i.e. R, 



i?in, the SR gain also 



measures the response of the non-linear system with re- 
spect to that of a linear system subject to the same de- 
terministic and stochastic forces. 

Additionally, in the absence of interaction between the 
bistable oscillators (the case 9 = 0), the collective SR 
gain G equals the SR gain of each independent oscillator 
Thus, by comparing the SR gain values of the collec- 
tive variable of a finite set of interacting oscillators with 
those observed in the case N = 1, we have a useful tool to 
highlight nonlinear effects that are a direct consequence 
of the coupling between the oscillators. 



III. FINITE SIZE DYNAMICS 

In this section, we define a set of stochastic processes, 
which we will refer to as fluctuating cumulants, in order 
to describe the dynamics of a finite system of coupled os- 
cillators in terms of a reduced number of variables. Then, 
by using Ito stochastic calculus, we derive the hierarchy 
of equations that these cumulants obey. A few approx- 
imation schemes are proposed for systems with a large 
but finite number of oscillators. Finally, we introduce a 
simple one- variable model in which the dynamics of S{t) 
is mimicked by using an effective potential. 

The fact that the infinite system {N = oo) is com- 
pletely deterministic, and the approximations described 
in this section are valid for large N, makes these methods 
especially appropriate to study the effect of fluctuations 
due to the flnite size of the system. 
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A. Fluctuating cumulants 



are 5— correlated noises with the following first moments 



Let us define the set of stochastic variables 

N 

N ■ 



1 ^ 



(11) 



i=l 



with k being a positive integer. We will refer to Mk as 
the fluctuating moment of order k. Note that Mi{t) = 0. 

In order to obtain a hierarchy of stochastic differen- 
tial equations for these variables, we need first to choose 
a convenient stochastic interpretation. The Langevin 
equations ([T]) are well defined and do not depend on the 
stochastic interpretation. Note, however, that the spe- 
cific form of the stochastic differential equations for the 
fluctuating moments do depend on the stochastic calcu- 
lus utilized. In the following, unless explicitly stated, Ito 
stochastic calculus is assumed. It is customary within 
this calculus to use a notation to express the stochastic 
differential equations in which there is no explicit men- 
tion to the white noises (see for example Q). In partic- 
ular, Eq. ID) would be written as 

9 ^ 

dx, ^[x^~xf + — J2i^J ~ 2;,) + F{t)]dt + (2Z3)i/2dB„ 

(12) 

where dBi , with i = 1, . . . ,N, is the differential of the 
Wiener process Bi{t) with properties 



dB,{t)dBj{t) = Sijdt. 



(13) 



The Gaussian white noise £,i{t) can be viewed as pro- 
portional to the derivative of Bi{t), — (2£>)^/^di?i/dt, 
though it is not an ordinary stochastic process but a gen- 
eralized process and requires a special formalism to be 
defined rigorously (see @ and references within). Here 
we will use both notations at convenience. 

Using Ito differentiation rules 9] we find the following 
stochastic differential equations for the fluctuating mo- 
ments 



k 



= (1-38^^ -9)Mk-Mk+2-3SMk+i 



+ (M3 + 3SM2)Mk-i + (fc - !)£> I 



N 



-r]Mk-i + fJ.k-1, 



(14) 



where 



1 ^ 



and 



f ^ 



(16) 



(17) 
(18) 

(19) 



{rj{t)) - (Aife(i))=0, 
2D 

mvit')) = ^sit'-t), 

2D 

{^ik{t)^ik'it')) - —6{t'-t){Mk+k'it)). 



Notice that the result is only obtained from ([T| 
when Ito calculus is assumed (see Appendix [XJ . Using 
Stratonovich calculus leads to a much more intricate ex- 
pression in which the approach proposed in this paper is 
not applicable. 

Additionally, note that (fT9|) implies that the processes 
/ifc are not uncorrelated. Rigorously, only r]{t) is a Gaus- 
sian process. Nevertheless, it can be shown that in the 
asymptotic limit of a very large number of oscillators, 

— > 00, all /Xfe tend towards a Gaussian behavior (see 
Appendix [B|. This property will allow us to rewrite 
Eq. (fH]) as a closed set of stochastic equations for the 
fluctuating variables in that limit. 

We can define a set of fluctuating cumulants by using 
the formula: 



n-l 
k=l 



-1)1 



k\{n - 1 - fc)! 



Kn-kMk. 



(20) 



Equation ([^D|) is the formula that relates the moments 
with the cumulant moments of a single- variable stochas- 
tic distribution. When the stochastic variable is Gaus- 
sian, all cumulants Kn, with n > 3, exactly vanish. A 
description in terms of cumulants is preferable because, 
unlike a descriptions with moments, it is expected that 
higher order cumulants are negligible in comparison with 
lower order cumulants, especially if the deviation with 
respect a Gaussian behavior is not very large. 

In terms of the fluctuating cumulants Kn, the first 
three equations of the hierarchy are 



S-S^- 3SK2 -K3+r] + F{t), 



— ^ = (1 - - 9)K2 - 35/^3 - 3Ki - Ki 



1 



S 

K_2 

2 



(21) 



-on — 

^ N 



K3 
3 



Ml, (22) 
(1 - 352 - 9)K3 - 3S{2Ki + Ki) 



-9K2K:!, - - K2ri + H2- 



(23) 



The noise terms 77 and /i^ vanish in the formal limit 
A^ — > 00, therefore leading to a deterministic hierarchy 
of equations for the fluctuating moments or cumulants. 
This deterministic hierarchy is equivalent to the non- 
linear hierarchy obtained by Desai and Zwanzig in [5|] for 
the cumulant moments of the process yi{t) = xx{t) — S{t). 
In contrast to the theory presented in Ref. which is 
based on the calculation of one-time expectation values, 
the fluctuating cumulant approach will allow us to study 
dynamical properties such as autocorrelation functions. 
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When TV is finite, the hierarchy of equations that Mk 
or K}. obeys is not closed, since the noise processes /j.^ de- 
pend on Mk in a non-trivial way. In the next subsections 
we present a few approximative schemes that overcome 
this difficulty for systems with a large number of coupled 
oscillators. 



B. Second order approximation 

If we retain the first two equations of (PT|l - (l^51) . neglect 
all fluctuating cumulants Kn with n > 3, and also neglect 
the term and the noise Hiit) in (122]) . we obtain a 
closed set of equations for the processes S{i) and K2{t): 

S = S-S^ -?,SK2 + r^ + F{t), 
2 



(1 - 35^ - e)K2 ~ 3Ki + D. 



(24) 



These set of equations was proposed by Pikovsky et al. in 
Ref. This truncation scheme has been called "the 
Gaussian approximation" because all fluctuating cumu- 
lants with order higher than the second one are neglected. 
There is no reason to expect a priori that these higher 
order cumulants can be neglected in any limit, other than 
the hope that their contribution is small. Note, in addi- 
tion, that in this scheme the (5-correlated noise is 
neglected without justification. 



C. Third order approximation 

Let us now focus in a third order truncation scheme. 
We will retain the three Eqs. (I^T1) - (P5)) . but consistently 
neglect K4 and K^. 

Since each iik{t) for k = 1,2,... is a Gaussian pro- 
cess in the lowest order in (see Appendix IB|1 . its 
probability distribution is completely determined by its 
first moments (|17p - (fT9|) . As mentioned before, the pro- 
cesses fik{t) SLie not independent of each other. Thus, it 
is preferable to express them in terms of a set of inde- 
pendent Gaussian noises rii{t) with zero mean and 



ID 

{m{i)m'(t'))^j^5u'5{t-t'), 

where l^V > and r]Q = rj. With the expansion 



(25) 



(26) 



we only need to select the coefficients Cki so that the cor- 
relations Uni) are satisfied. This can be achieved by using 
the Gram-Schmidt ortho-normalization method. The re- 
sult for the first two terms is 



Ml 



{K2)^rii 
{K2)V- 



(27) 



\{K4){K2)+2{K2)'-{K-i)^\^ 



m- (28) 



Note, however, that in these expressions the coefficients 
Cfei appear as functions of the average values of the fluc- 
tuating cumulants Kn- Thus, if we plan to solve the Eqs. 
(l2T])-(l23l) using (I27l)-(l28l), wc would have to consider the 
equation of motion for (Kn) 5] and solve the whole set 
of equations self-consistently. Alternatively, we could use 
a slightly different version of (P7)) - ([^51) in which the ex- 
pected values (Kn) are replaced by Kn. This way, the 
correlations for the first noise terms are also iden- 
tically satisfied, and the fact that the fluctuating cumu- 
lants become deterministic in the limit N —>■ 00 guaran- 
tees that the proposed expressions for the noises /ife(t) 
are Gaussian in the lowest order in N~^. Since all Gaus- 
sian processes are completely determined by its first two 
moments, both methods to generate the noise sources /ife 
are mathematically equivalent in the asymptotic limit of 
large N, though the later is more physically appealing 
because in this case the instantaneous value of the noise 
source /ife(i) in one trajectory does not depend on aver- 
ages over trajectories but on single-trajectory values. 

Using this later procedure, the following stochastic dif- 
ferential equations with multiplicative noise are obtained 



S 

K2 
2 

Ks 
3 



= S -S^ - 3SK2 - + + F{t), 



(1 - 352 - e)K2 - 3SK3 
+D + \K2\irji, 



iK^ 



(1 - 352 



)K3 - &SK'^ - 9K2K3 



K:,7^i + \2Kl~Kl\2 7j2 
IKol^ 



(29) 



■m 



This method has also the advantage that the system of 
equations (j29p can be solved numerically using standard 
stochastic algorithms [l^. It represents a third order 
approximation scheme. Notice that this scheme can be 
applied in a straightforward way to obtain the corre- 
sponding equations of an arbitrary order truncation of 
the fluctuating cumulant hierarchy. 



D. Effective potential 

As we increase the truncation order of the fluctuating 
cumulant hierarchy, as we have discussed above, a more 
accurate approximation is obtained. However, the num- 
ber of equations is also increased. On the other hand, 
we may wonder how good a description based on a single 
differential equation is. The aim is to derive a simplified 
model that may not mimic quantitatively but qualita- 
tively the coupled system dynamics, in addition to being 
more amenable to analytical treatment. 

Here we consider the following single stochastic equa- 
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tion 

S=-U:,f,{S)+v + F{t), (30) 

where ri{t) is the Gaussian white noise defined by l|17p- 
(fTS]) . and Ucs{S) an effective potential to be specified. 

We can determine the effective potential uniquely by 
requiring the model to reproduce the equilibrium prop- 
erties of the original system. The stationary probability 
density Peq{S) of the Langevin Eq. ([501) i'^ the absence 
of external driving is given by [Tlj 

Pcq{S) = Z exp I 1, (31) 

where Z is a normalization constant. Therefore, by in- 
verting (j31[) . we find an expression for the effective po- 
tential up to an additive constant c, 

C^cff(^) = -^lnPeq(^)+c. (32) 

In Ref. 0, Desai and Zwanzig presented an analytical 
expression for the equilibrium density Pcq(»S') by retaining 
the leading term in the asymptotic expansion for large N. 
We will refer as U^^\s) to the corresponding effective 
potential. This analytical solution shows that for a given 
value of d there exists a. D = Dc such that for values 
D greater than D^, the effective potential uj,^\s) is 
monostable with a minimum located at 5*0 = 0. For 
< D < Dc, the effective potential is bistable with two 
minima at iSo, being 5*0 a function of 9 and D. Figure [1] 
depicts this situation for a system with 9 = 0.5. The 
calculated critical noise for this value of 9 is Dc ^ 0.2645 

For a system with a finite size N, we can calculate nu- 
merically Pcq{S) by simulating the Langevin equations 
(fT|) and computing the histogram of the collective vari- 
able S{t) after a sufficiently large time when the system 
has equilibrated. Figure [T] shows the resulting effective 
potential U^g\s) for a system with iV = 10 oscillators. 
It can be seen that the deviations with respect to the 
infinite size potential U^g\S) are very small, even for 
such a small system. 



IV. STOCHASTIC RESONANCE REVISITED 

In this section, we compare numerically the predictions 
of the effective models presented in Sec. Illll in the frame- 
work of SR. The simplified character of these models will 
allow us to explain in intuitive terms the highly nonlinear 
effects observed in the stochastic resonance quantifiers 

We will restrict our study to a periodic rectangular 
driving force, 

F{t) = (-1)"(*U, (33) 




FIG. 1: Effective potential for the simple Langevin model 
(|30p for two systems with 6 — 0.5. Top panel corresponds to 
D = 0.08 (< Dc ~ 0.2645) and the bottom panel to D = 0.4 
(> Dc)- The solid lines depict the analytical solution given in 
Ref. [a], whereas the dotted lines correspond to the effective 
potential obtained using the simulation method described in 




FIG. 2: Effective dynamics for a coupled system with A*' = 
10 oscillators subject to a rectangular driving of amplitude 
A = 0.3 and frequency fl — 0.01. Several stochastic reso- 
nance quantifiers are depicted vs the noise strength D. From 
top- left to right-bottom: the numerator of the SNR (Qu), 
the denominator (Qi), the output SNR (-Rout) and the gain 
(G). The solid lines depict the numerical solution of the full 
Langevin dynamics ([l]). The dotted and dashed lines corre- 
spond to the fluctuating cumulants approach truncated at the 
second (|24|) and third order (|29|) . respectively. The results for 
the effective potential approach (|30p are depicted by triangles 
pointing upwards (U^^^) and downwards ((7^^"'). 

where n{t) — [2t/T\, [zj being the floor function of z. 
The input SNR for forces of this type can be readily cal- 
culated as i?in ~ 4:A'^N/{ttD). In all cases reported here 
the coupling strength is fixed to 9 = 0.5, the driving fre- 
quency Q = 2tt/T to Q — 0.01, and the driving amplitude 
to ^ = 0.3. This amplitude is subthreshold in the sense 
that the driving force (|33p cannot induce sustained oscil- 
lations between the different attractors of the dynamics 



in the absence of noise (i.e., for D = 0). 

The stochastic differential equations presented in the 
preceding sections were solved numerically by using weak 
predictor-corrector algorithms of order 2.0 [lO|. 

Figure [2] shows several SR quantifiers as a function of 
the noise strength D for a coupled system with = 10 
oscillators. A strong amplification of the collective re- 
sponse is observed, with SR gains reaching very large val- 
ues, especially when compared with uncoupled systems 
subject to the same input signals (see Ref. [13 )■ These 
findings were first reported in Since the numerator 
Qu of the SNR remains of the same order of magnitude 
for the range of noise strength values plotted, the large 
values of the SR gain are mainly due to the reduction 
of a few orders of magnitude of the denominator Qi, as 
shown in the top-right panel of Fig. [21 

The SR quantifiers obtained with the effective poten- 
tial model described by Eq. (fSO]) are depicted by triangles 
in Fig. [21 Triangles pointing upwards correspond to the 
effective potential C^'^^ in the asymptotic limit N — )■ oo, 
whereas triangles pointing downwards correspond to the 
effective potential uj,^^ computed numerically for a sys- 
tem with N = 10. In Fig. [21 it can be seen that the 
later leads to a better agreement than the former for 
Qu due to the small but appreciable discrepancies ob- 
served in Fig. [1] However, no significant improvement is 
seen in the rest of the quantifiers: Qi, the SNR, and the 
gain. In general, the effective potential approach is able 
to describe qualitatively the phenomenon, displaying a 
non-monotonic behavior with a maximum at about the 
same value of the noise strength D than the original sys- 
tem. Nevertheless, quantitatively the agreement is not so 
good, showing a consistent underestimation of the noise 
term Qi by roughly a factor of 2. The fact that this ap- 
proach underestimates the fiuctuations of the collective 
variable is easy to understand if one takes into account 
that the effective potential is a mean-field-like idealiza- 
tion in which the real discrete, more noisy, interaction is 
replaced by a smoothed potential. 

A slightly better quantitative agreement is obtained 
with the Gaussian approximation described by Eq. (j24p . 
i.e. the fiuctuating cumulants approach truncated at 
the second order, which is depicted in Fig. [21 by dotted 
lines. It can be seen that the SNR is in better agree- 
ment, though the SR gain around the maximum has not 
been improved significantly overall. A considerably en- 
hanced agreement is achieved by the third order approach 
in Eq. (I29p . which is represented in Fig.[21by dashed lines. 
It can be seen that this third order approximation slightly 
underestimates the noise term Qi for large enough val- 
ues of D. This is what one would expect, because this 
method neglects the higher order cumulants K4 and A'5 
in Eq. ^ [see Eqs. ((221) "(123)], which, if present, would 
increase the fluctuations of the lower order cumulants. 

Figure [3l confirms the above discussed behavior of the 
simplified models for a larger system with A^ = 30. 
Again, the effective potential theory and the Gaussian 
approximation provide mainly a qualitative picture, with 




FIG. 3: Same as in Fig. [2] but for a system with iV = 30. 



quantitative predictions within the same order of magni- 
tude. The best quantitative agreement is also observed to 
be given by the third-order fluctuating cumulant scheme. 
The main differences with respect to the smaller system 
discussed before are in the quantity Qu, which is pro- 
portional to the spectral amplification. In this case the 
system is large enough so that very small differences are 
observed between the effective potential uj,^^'' and uj,^\ 
They both provide data in very good agreement with 
the original system data. In addition, notice that the 
Gaussian approximation data for Qu deviates apprecia- 
bly from the system data for large enough values of D. 
The fact that the third order approximation leads to a 
good agreement for Qu indicates that the third cumulant 
plays an important role for the spectral amplification at 
these noise strength values. 

Finally, we now use these simplified models to explain 
the very large gain values observed in globally coupled 
bistable systems 0, and particularly those observed 
in the bottom-right panels of Figs. [21 and [31 These gain 
values are much larger than those observed in uncoupled 
or isolated bistable systems subject to the same rectan- 
gular input signals [12], and thus are due to the inter- 
action between the oscillators. For the sake of simplic- 
ity we will use the effective potential approximation in 
the following discussion, though similar arguments can 
be used within the other schemes presented in Sec. IIIII 
The only noisy term in Eq. (j30p is ri{t), which has a 
strength of D/N. Thus, we can fix the individual noise 
strength D and still being able to get rid of the noise 
by considering the limit N —>■ 00. This way Eq. (|30p 
becomes deterministic and we can apply the concept of 
static threshold for a finite value oi D. A simple anal- 
ysis of the infinite size potential shows that a constant 
driving of ^ — 0.3 is able to remove one of the two at- 
tractors of the dynamics for systems with noise strength 
values D larger than Dc{A = 0.3) « 0.02 (and, of course, 
smaller than Dc{A = Q) ^ Dc~ 0.2645, the noise value 
where the effective potential turns monostable in the ab- 
sence of driving). Therefore, most of the data points 
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FIG. 4: Subthreshold dynamics. Same stochastic resonance 
quantifiers as in Fig. [2] but as a function of the system size 
for a fixed noise strength D — 0.017. Squares depict the 
simulation data corresponding to the full dynamics ([l]). Lines 
are a suide to the eye. 




' ff. I . I . I . I . I . I . I .1 - A ,1,1,1,1,1,1,1,' 
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FIG. 5: Suprathreshold dynamics. Same as in Fig. |4] but for 
D = 0.2. 

in Figs. [2] and [3] correspond to suprathreshold dynamics 
when viewed from the perspective of the Langevin equa- 
tion ([50)1 . The response is expected to be more amplified 
with suprathreshold signals than with subthreshold sig- 
nals, because, for the first ones, the presence of noise is 
not necessary in order to produce jumps between the two 
locations of the time-dependent attractor. For instance, 
only gains larger than unity has been found for an iso- 
lated bistable system subject to monochromatic signals 
when the driving amplitude is suprathreshold p^ . with 
gain values reaching a few tenths above unity. In fact, the 
above consideration may well explain why gains larger 
than unity (also a few tenths above unity) are found for 
globally coupled bistable systems subject to monochro- 
matic signals the collective variable dynamics is ef- 
fectively suprathreshold. 

A numerical analysis of the deterministic version (i.e. 
without the noise terms rjk) of Eqs. (|29p governing the 



third order approximation confirms that under a con- 
stant driving of A — 0.3, there is a transition at about 
0.02 between a situation in which the system presents 
two attractors (subthreshold dynamics) and only one 
(suprathreshold dynamics). In fact, this analysis can also 
be carried out with the fluctuating cumulant theory pre- 
sented in Sec. IIIII with an arbitrary order of truncation, 
and thus, with an arbitrary order of accuracy. 

Let us illustrate the above discussion by considering 
two systems: one with a noise strength value D = 0.017, 
which is just below the transition value Dc{A = 0.3); and 
one with D = 0.2, well above that transition value. Fig- 
ures m and O shows the behavior of the stochastic reso- 
nance quantifiers as a function of the system size N for 
these two systems under the same rectangular driving 
with A = 0.3 and fl — 0.01. It can be seen that, for 
large enough A^, the size of the fluctuations (Qi) is al- 
ways reduced as N is increased, as expected. However, 
both systems display a very different behavior. In Fig. [4] 
the SR quantifiers Qu, SNR and the gain are observed to 
decrease monotonically as N is increased. This is what is 
expected for a system under subthreshold dynamics, be- 
cause the driving force alone is unable to provoke jumps 
between the attractors and needs the presence of fluc- 
tuations. In the limit N —> oo, the fluctuations, and 
thus the jumps, are completely suppressed. On the other 
hand, a very different behavior is observed in Fig. [5] for 
a noise strength value above the transition value corre- 
sponding to this driving amplitude. In this case, there is 
only one attractor, which is displaced to the positive or 
negative axis around the origin according to the instanta- 
neous sign of the driving force. Here the presence of finite 
size fluctuations is only a nuisance to the driven oscilla- 
tions of the collective variable. As a result, the spectral 
amplification and the SR gain grow with N toward the 
finite values which correspond to the infinite system. 

Finally, let us notice that the phenomenon of system 
size resonance 0], i.e. the non-monotonic behavior of 
the SR quantifiers when plotted as a function of N, 
could only appear in a parameter region in which the 
system dynamics is subthreshold. Nevertheless, this phe- 
nomenon is not observed in Fig.|3) This is due to the fact 
that the strength of the fluctuations is very small even 
for small values of the system size N, already below the 
optimal noise value in which the response is maximized. 
Any further increase of the system size reduces the fluc- 
tuations and, thus, the jumps between the attractors. 



V. CONCLUSIONS 

Starting from the Langevin equations defining the 
model system, we have derived, using Ito stochastic cal- 
culus, a hierarchy of exact stochastic differential equa- 
tions for a set of fiuctuating cumulant variables, defined 
by using the arithmetic mean over all oscillators. Due 
to the useful mathematical properties of Ito stochastic 
calculus, the hierarchy contains noise terms with simple 
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autocorrelation properties. Furthermore, the approach 
proposed in this paper for the fluctuating cumulants is 
not directly applicable with Stratonovich calculus. In the 
limit of an infinite number of oscillators, the whole hierar- 
chy reduces to the one obtained by Desai and Zwanzig [5] 
for the expected values of the cumulants. In contrast to 
the theory presented in Ref. [5| , the fluctuating cumulant 
approach allows us to study a wide kind of collective dy- 
namical properties like autocorrelation functions or the 
SNR, in addition to effects due to finite size fluctuations. 

Nevertheless, the noise terms that appear in the exact 
hierarchy for the fluctuating cumulant variables depend 
in a complicated way on the fluctuating cumulants and 
approximations have to be taken in order to obtain a 
closed set of stochastic differential equations. Here it is 
shown that this difficulty can be overcome in the asymp- 
totic limit of a very large number of oscillators. How- 
ever, one still has the inconvenience of dealing with a 
hierarchy with an infinite number of equations, and a 
truncation scheme is desirable. A Gaussian approxima- 
tion was proposed by Pikovsky et al. in Ref. [1], and 
here it is presented as a second-order truncation scheme 
of the fluctuating cumulant hierarchy. In addition, an 
arbritrary-order truncation method is proposed, with ex- 
plicit expressions given for the third order only. This 
third order approach turns out to provide the best quan- 
titative agreement with the SR data. Finally, a rather 
simple approach based on a single variable and the use 
of an effective potential is proposed. 

The spectral amplification of the collective variable as 
a function of the noise strength D of systems with = 10 
and A'' = 30 bistable oscillators is found to be in good 
agreement with the predictions given by all the approxi- 
mations, though small but appreciable systematic devia- 
tions are observed for the Gaussian approximation for a 
system with A = 30 oscillators. Additionally, the effec- 
tive potential theory and the Gaussian approximation do 
not account well for the SNR or the SR gain of the col- 
lective variable, though the data is within the same order 
of magnitude. A systematic underestimation of the fluc- 
tuations of the collective variable is done by the effective 
potential approach due to the mean-field like character of 
this simplified theory. The best quantitative agreement 
of the SNR and the SR gain is given by the third-order 
fluctuating cumulant theory, although a small systematic 
underestimation of the fluctuations is still observed with 
this third order theory due to the neglect of higher order 
fluctuating cumulants. 

Furthermore, using any of the approximations pre- 
sented, we are able to explain the very large gain val- 
ues observed in Refs. Specifically, it is shown that 
the driving amplitudes used are suprathreshold from the 
point of view of the effective dynamics in the range of 
noise strength values utilized in those works, i.e., there is 
only one attractor of the dynamics under the presence of 
the driving force, and this attractor oscillates following 
the driving force. Simulation results, showing several SR 
quantifiers as a function of the system size N , confirm 



this behavior. 

This situation resembles the effect of a high-frequency 
signal in an isolated bistable system. In the later case, 
the high-frequency signal can be removed from the de- 
scription by means of an effective bistable potential 
with modified parameters, with the consequence that 
previously subthreshold driving amplitudes can become 
suprathreshold from the point of view of the effective po- 
tential In contrast, the effective dynamics induced 
by the high-frequency signal has been shown to provoke 
the opposite effect on an excitable system, beiiig able to 
suppress the excitable character of the system [15| . This 
suggests that much work is needed in order to extend the 
present analysis to finite sets of coupled excitable systems 
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APPENDIX A: DERIVATION OF EQ. (fT9l) 

In the notation commonly utilized within the frame- 
work of Ito calculus, Eq. can be expressed as 

2D 

dTk{t)dVi{t) = —Mk+i{t)dt, (Al) 



ATk{t)dTi{t') =0 for<^<'. 



(A2) 



where Tk{t) = /J dT/ifc(T). To prove these equations we 
start from the definition (fT51) to arrive at 



dTk{t) = ^[y.(t)]Mi?,(i), (A3) 



A^ 



i=i 



where yi{t) = a;,(i) - S{t). Thus 
2D 



dTk{t) ATi{t) 



A2 

2D 

A2 

2D 



i 

Mk+i{t)dt, (A4) 



where we have used p^ . Similarly, using the fact that 
the Wiener processes Bi{t) have independent increments, 
i.e. dBj{t) dBj{t') = for t ^ t' , and that the increments 
dBj{t') are independent of yi{t) at times t' > t, Eq. (|A2|) 
is readily proven. 
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In order to clarify the advantages of Ito calculus in 
the context of this paper, let us now compute the auto- 
correlation of Hk{t) by using Stratonovich calculus. In 
this case we are entitled to utilize the usual rules of dif- 
ferentiation of deterministic calculus. To that aim, the 
Novikov-Furutsu theorem [l7|, d, d states that if ^(t) 
is a Gaussian white noise with zero mean and autocor- 
relation {^{t)£,{s)) = 2D5{t — s), then for any functional 
g[^] we have 



mm) = 



Mmas))M) = 2Dm)^ (A5) 



'my 



where 5g[^]/5£,{t) denotes the functional derivative of 
g[£]. Thus, assuming t <t', 



2D 



m)yAt 



sm') 
-) + 



(A6) 



Using repeatedly the Novikov-Furutsu theorem and the 
fact that 5yj{t)/6^i{t) = {l/2)Sij, we arrive at Eq. (HH) 
plus the following two extra terms on the righthand side 
of Eq. (fT9| 



D^ik{Mkm)Mim)) 



2DH 



sm 



. (A7) 

These extra terms make the problem much more difficult 
to deal with. 



third moment of /ife(t). 



{^lk{tl)^m■2)^^k{h)) = o. 



(Bl) 



and all odd moments of /ifc {t) vanish. If /ife it) were Gaus- 
sian, all cumulants higher than the second should be zero. 
This requires all odd moments of /ife(t) to vanish but also 
a specific functional form of the even moments flj]j. For 
instance, were /ife(i) a Gaussian process, the fourth mo- 
ment {iJ.k{ti)lJ-k{t2)iJ'k{h)lJ.k{ti)) would be given by 



{tl)^J■k {h)){^.k {h)^.k {U)) 
+ {nk{ti)^.k{h)){^ik{t2)^'k{ti)) 

+ {iXk{tl)iik{ti)){lJLk{t2)ilk{t-i))- (B2) 

Instead, a simple calculation shows that for a finite N 
the fourth moment is 



{lJ-k{ti)lJ-k{t2)i^k{h)lJ-k{t'i)) 



2D\- 



N J 

{M2k{tl)M2k{t2))5{tl - U)5{t2 - h) 
f {M2k{tl)M2k{t2))&{tl - h)6{t2 - U) 

f {M2k{t2)M2k{h))5{ti - t2)5{U - U)] . (B3) 



Therefore, by considering it is clear that (|B2[) does 
not equal (|B3|) unless 



{M2k{tl)M2kit2)) = {M2kitl)){M2kit2)). 



(B4) 



APPENDIX B: GAUSSIAN NOISES IN THE 
LIMIT OF LARGE NUMBER OF OSCILLATORS 

In this appendix we show that the process /ifc (i) tends 
to a Gaussian behavior as TV — *■ c>o. First note that the 



This identity is asymptotically correct in the limit N — > 
oo, because then all fluctuating moments become deter- 
ministic. For a large enough TV, Eq. (|B4|) will hold as 
a good approximation, showing a Gaussian behavior in 
the lowest order of a expansion. Clearly, similar 

considerations apply to other even moments of /ife(i). 
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